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ABSTRACT:  We  studied  the  thermal  decomposition  and  subsequent 
reaction  of  the  energetic  material  nitromethane  (CH3N02)  using  molec¬ 
ular  dynamics  with  ReaxFF,  a  first  principles-based  reactive  force  field.  We 
characterize  the  chemistry  of  liquid  and  solid  nitromethane  at  high 
temperatures  (2000—3000  K)  and  density  1.97  g/cm3  for  times  up  to 
200  ps.  At  T  =  3000  K  the  first  reaction  in  the  decomposition  of 
nitromethane  is  an  intermolecular  proton  transfer  leading  to  CH3NOOH 
and  CH2N02.  For  lower  temperatures  (T  =  2500  and  2000  K)  the  first 
reaction  during  decomposition  is  often  an  isomerization  reaction  involving  the  scission  of  the  C— N  bond  the  formation  of  a  C— O 
bond  to  form  methyl  nitrate  (CH3ONO).  Also  at  very  early  times  we  observe  intramolecular  proton  transfer  events.  The  main 
product  of  these  reactions  is  H20  which  starts  forming  following  those  initiation  steps.  The  appearance  of  HzO  marks  the  beginning 
of  the  exothermic  chemistry.  Recent  quantum-mechanics-based  molecular  dynamics  simulations  on  the  chemical  reactions  and  time 
scales  for  decomposition  of  a  crystalline  sample  heated  to  T  =  3000  K  for  a  few  picoseconds  are  in  excellent  agreement  with  our 
results,  providing  an  important,  direct  validation  of  ReaxFF. 


■  INTRODUCTION 

A  molecular-level  understanding  of  condensed-phase  chem¬ 
istry  is  among  the  central  challenges  of  the  atomistic  modeling 
community;  such  a  fundamental  understanding  is  critical  in  many 
areas  including  materials  for  the  energy  sector  and  high-energy 
density  materials.  Ab  initio  quantum  mechanical  (QM)  methods 
to  obtain  first  principles  predictions  of  the  mechanisms, 
barriers,  and  rates  of  unimolecular  (see  for  example  ref  l)  and 
bimolecular2  reactions  have  dramatically  improved  our  under¬ 
standing  of  chemical  reactions  relevant  to  many  fields  including 
energetic  materials1,2  and  catalysis."’  Such  QM  studies  are 
computationally  intensive  and  are  typically  limited  to  relatively 
small  systems  and  short  times  (roughly  ~100  atoms  for  a  few 
picoseconds).  Furthermore,  QM  studies  usually  focus  on  gas 
phase  phenomena;  thus,  pressure  effects  and  multimolecular 
reactions  can  only  be  taken  into  account  approximately. 
An  exception  are  recent  QM  based  molecular  dynamics  (MD) 
studies  of  high  temperature  decomposition  of  energetic  materials4-6 
that  follow  condensed-phase  dynamics  of  the  system  in  full  atomistic 
detail.  The  computational  cost  of  QM  limits  most  simulations  to  small 
systems  over  the  course  of  a  few  picoseconds.  A  recently  developed 
method  allows  coupling  of  a  small  QM  region  with  a  continuum 
description  of  a  much  larger  system.7,8  This  illuminates  some  aspects 
of  the  physics  and  chemistry,  but  full  scale  molecular  dynamics  on 
larger  systems  are  still  needed. 

To  achieve  longer  simulation  times  in  larger  systems,  it  is 
possible  to  replace  the  QM  description  with  a  force  field  (FF) 
that  describes  atomic  interactions  in  a  computationally  efficient 
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manner.  This  has  allowed  MD  simulations  for  up  to  microse¬ 
conds  and  is  routinely  used  for  tens  of  nanoseconds.  A  severe 
limitation  of  this  approach  has  been  that  such  classical  FFs  have 
been  incapable  of  describing  correctly  the  breaking  and  forming 
of  covalent  bonds,  that  is,  reaction  chemistry.  Over  the  past  few 
years  we  have  been  developing  a  new  generation  of  FF  (denoted 
ReaxFF)  aimed  at  accurately  describing,  in  a  computationally 
efficient  manner,  the  chemistry  of  a  wide  range  of  materials. 
ReaxFF  uses  general  relationships  to  describe  how  charge  flows 
between  atoms  and  how  the  valence  energies  depend  on 
geometry.  Thus,  the  valence  energies  are  described  in  terms  of 
bond  orders  that  are  uniquely  determined  from  the  atomic 
positions  (no  predetermined  assignment  of  bond  types),  and 
the  parameters  used  in  ReaxFF  are  obtained  from  QM  calcula¬ 
tions  on  structures  and  reactions  for  a  large  set  of  systems. 
ReaxFF  has  now  been  used  successfully  for  a  number  of 
allowed  and  forbidden  reactions  of  hydrocarbons,9  oxidation  of 
metals  such  as  Si,  Al,  and  Pt, 10,11  and  for  shock12and  thermal1"’ 
induced  decomposition  the  energetic  material  RDX  [cyclic 
(CH2N(N02))3J.  The  use  of  this  type  of  force  field  to  predict 
complex  processes  in  real  materials  (often  under  high  tempera¬ 
tures  and  pressure)  is  growing  at  a  rapid  pace,  and  there  is  a 
serious  need  for  a  critical  assessment  of  their  accuracy.  Among 
the  many  possible  applications  of  ReaxFF,  the  decomposition 
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and  reaction  of  energetic  materials  are  among  the  most  challen¬ 
ging  ones  due  to  the  complex  chemistry  involved  and  the  extreme 
conditions  of  temperature  and  pressure  involved. 

In  this  paper  we  use  ReaxFF  with  MD  to  study  the  thermal- 
induced  decomposition  and  subsequent  reactions  of  the  ener¬ 
getic  material  nitromethane  (CH3N02)  at  high  temperatures 
(from  2000  to  3000  K)  and  high  density  p  =  1.97  (normal  density 
is  1.13  g/cm3)  for  times  up  to  200  ps.  We  report  here  the  thermal 
decomposition  of  liquid  and  crystalline  samples  of  nitromethane 
that  are  instantaneously  heated  to  the  desired  final  temperature. 
Reactive  MD  simulations  provide  a  molecular-level  description 
of  the  chemistry  and  mechanical  processes  in  energetic  materi¬ 
als  making  no  assumptions  or  simplifications  other  than  those 
intrinsic  to  the  force  field  description  of  the  interactions,  the 
fact  that  we  use  classical  mechanics  to  describe  atomic 
dynamics,  and  the  size  of  the  simulation  cell.  Our  simulations 
show  that  the  first  chemical  events  involve  an  intramolecular 
proton  transfer  leading  to  CH2N02  and  CH3NOOH  or  an 
isomerization  reaction  leading  to  CH3ONO.  This  is  the  first 
time  that  such  isomerization  has  been  observed  in  condensed- 
phase  simulations.  We  also  find  intramolecular  proton  transfer 
at  early  times  (within  about  a  picosecond  for  T  =  3000  K).  After 
these  initial  events,  N— O  bond  breaking  leads  to  the  formation 
of  H20  and  the  initiation  of  the  exothermic  chemistry.  At  this 
point  the  potential  energy  decreases  dramatically  as  the 
pressure  increases  while  the  average  molecular  weight  of  the 
system  decreases.  We  choose  nitromethane  due  to  the  exten¬ 
sive  experimental  characterization  of  this  material  and  also  to 
compare  the  ReaxFF  results  with  recent  work  by  Manaa  and 
collaborators  who  used  QM-based  MD  to  heat  a  crystal  to  3000 
K  for  time  scales  of  around  1  ps.  We  find  very  good  agreement 
between  our  ReaxFF  simulations  and  QM-based  ones  both 
regarding  nature  of  the  initial  chemical  reactions  as  well  as  their 
time  scales;  this  provides  an  important  validation  of  the 
ReaxFF  description  of  reactive  atomic  interactions  of  CHNO 
materials. 

After  a  brief  description  of  relevant  previous  work  (section  II), 
in  section  III  we  describe  ReaxFF  and  our  MD  simulations.  In 
section  IV  we  present  the  energetics  of  decomposition  and 
pressures  involved.  Section  V  describes  the  key  chemical  events, 
and  section  VI  focuses  on  the  atomic  diffusivity  of  reactants  and 
products.  Finally,  section  VII  draws  conclusions. 

■  BACKGROUND 

Nitromethane  is  one  of  the  most  extensively  studied  energetic 
materials  both  experimentally  and  theoretically.  It  is  generally 
agreed  that  under  low  pressure  conditions  the  first  reaction  in  the 
decomposition  of  NM  is  the  homolytic  breaking  of  the  C— N 
bond;  this  reaction,  having  a  positive  volume  of  activation,  does 
not  play  an  important  role  under  high  pressures. 

Several  observations  (see  for  example  ref  14)  indicate  that  the 
presence  of  the  aci  ion  form  of  NM  (CH2N02)  has  an  important 
effect  in  accelerating  the  chemical  reactions.  Engelke  and  colla¬ 
borators  showed  that  the  aci  ion  form  is  the  only  new  species  in 
base-sensitized  NM.15  Recent  infrared  spectroscopy  experiments 
show  that  intramolecular  proton  transfer  leading  to  the  aci  form 
is  among  the  first  events  in  the  room  temperature  decomposition 
of  NM  at  high  pressures.16  In  addition,  the  decrease  of  reaction 
time  with  increasing  pressure17  has  also  been  attributed  to  the 
increased  population  of  the  aci  ion  with  pressure.18  Time-of-flight 
mass  spectroscopy  studies  of  the  reaction  zone  of  detonating 


base-sensitized  NM  show  species  resulting  from  condensation 
reactions  involving  two  or  more  NM  molecules. 

From  a  theoretical  point  of  view  the  recent  development  of 
partial  bond-order-based  reactive  force  fields  opens  the  oppor¬ 
tunity  to  characterize  the  chemistry  of  energetic  materials  under 
conditions  approaching  those  in  experiments. 

Developing  a  fundamental  understanding  of  energetic  materi¬ 
als  requires  characterizing  complex,  inter-related  mechanical  and 
chemical  processes  at  extreme  conditions  of  temperature  and 
pressure  and  is  essential  for  the  development  of  physics-based, 
predictive  models  of  energetic  materials  that  can  aid  the  design 
of  safer  materials  with  improved  properties.  A  molecular-level 
characterization  of  the  processes  of  decomposition  and  reaction 
of  condensed-phase  energetic  materials  is  a  very  challenging 
problem,  and  force  fields  for  MD  simulations  capable  of  describ¬ 
ing  chemical  reactions9,19  are  playing  an  increasingly  important 
role.  Large-scale  MD  with  bond-order-dependent  reactive  po¬ 
tentials  that  exhibit  simple  exothermic  chemistry20  have  provided 
very  valuable  information  about  generic  properties  of  energetic 
materials  such  as  the  equilibrium  structure  of  steady  state 
detonation  waves,  '  role  of  defects,  '  the  existence  of  a  failure 
diameter,24  and  desensitization.25  The  recent  development  of 
reactive  potentials9,19  that  can  describe  the  chemistry  of  real 
materials  (including  HE)  has  enabled  the  characterization  of 
shock-induced  chemistry12  and  thermal  decomposition  in  a  wide 
temperature  and  pressure  range.12  However,  before  quantitative 
conclusions  on  real  materials  can  be  drawn  from  these  simulations 
a  critical  assessment  of  their  accuracy  needs  to  be  performed. 

■  COMPUTATIONAL  METHODS 

Reactive  Force  Field  ReaxFF.  The  key  concepts  that  enable 
ReaxFF  to  describe  complex  chemical  reactions  are  the  use  of  ( 1 ) 
environment  dependent  atomic  charges  (updated  at  every  step 
during  the  dynamics)  to  describe  electrostatics,  (2)  partial  bond 
orders  (a  many  body  contribution  based  purely  on  atomic 
positions)  to  describe  covalent  interactions,  and  (3)  two-body 
terms  (between  all  atoms)  to  account  for  London  dispersion 
attraction  and  short-range  Pauli  repulsion  (often  denoted  as  van 
der  Waals  terms). 

We  have  recently  described  the  energy  expression  of  ReaxFF 
and  the  training  set  used  to  parametrize  the  potential  for  RDX 
simulations.  This  FF  has  been  recently  improved  with  the 
addition  of  QM  results  on  chemical  reactions  involving  nitro¬ 
methane  and  vibrational  information  (frequencies  and  modes)  of 
small  molecules  to  the  training  set.  The  details  of  the  new 
chemical  reactions  used  for  the  ReaxFF  parametrization  were 
published  in  ref  26. 

Molecular  Dynamics  Simulations.  We  use  ReaxFF  with  MD 
to  simulate  the  decomposition  of  liquid  and  crystalline  samples  of 
nitromethane  at  temperatures  from  T  =  2000  K  to  T  =  3000  K 
and  density  p  =  1.97  g/cm3  (corresponding  to  a  volume  75%  of 
the  zero  pressure  volume  of  the  crystal).  The  crystalline  samples 
contain  16  molecules  and  were  obtained  by  replicating  the  unit 
cell  two  times  in  the  y  and  z  directions  followed  by  a  thermaliza- 
tion  at  T  =  10  K.  This  same  structure  was  used  in  ref  4  with  QM- 
based  MD  to  characterize  the  initial  chemistry  of  nitromethane  at 
T  =  3000  K.  The  liquid  samples  contain  32  nitromethane 
molecules  and  are  built  using  the  nonreactive  force  field 
DREIDING;27  we  perform  high  temperature  (T  =  500  K)  MD 
simulations  starting  with  a  low  density  simulation  cell  (where  the 
molecules  have  very  high  mobility)  and  compress  it  until  the 
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Figure  1.  Time  evolution  of  the  potential  energy  (in  keal/mol  per  nitromethane  molecule)  for  temperatures  3000  K  (a),  2500  K  (b),  and  2000  K  (c). 
Triangles  show  the  time  corresponding  to  the  first  N— -O  bond  breaking  event. 


desired  density  (1.97  g/cm3)  is  achieved.  This  structure  is  then 
quenched  with  isothermal  isochoric  MD  using  the  ReaxFF 
potential  at  T  =  10  K  to  obtain  the  first  liquid  sample  (denoted 
liquid  l).  We  then  run  40  ps  NVT  MD  at  T  =  500  K  with  ReaxFF 
and  rethermalize  the  resulting  structure  at  T  =  10  K  to  obtain  a 
second,  uncorrelated,  liquid  structure  (liquid  2). 

We  simulate  the  decomposition  of  all  three  structures  using 
isothermal  isochoric  MD  (NVT  ensemble)  using  a  Berendsen 
thermostat  with  a  coupling  constant  of  200  fs.  The  initial 
velocities  are  chosen  from  the  Maxwell— Boltzmann  distribution 
with  temperature  two  times  the  desired  temperature.  Since  the 
initial  potential  energy  is  close  to  the  ground  state,  the  kinetic 
energy  decreases  to  roughly  half  its  initial  value  very  quickly  (few 
tens  of  femtoseconds).  We  also  perform  an  annealing  simulation 
of  a  crystalline  sample  at  p  =  2.20  g/ cm3.  This  sample  was  heated 
from  T  =  1000  to  4000  K  to  compare  the  resulting  pressure 
profile  with  QM  results  reported  in  ref  4. 


Figure  2.  Time  evolution  of  the  pressure  for  the  liquid  structure  at  three 
temperatures:  T  =  3000  K,  T  =  2500  K,  and  T  =  2500  K. 


■  INTERNAL  ENERGY  AND  PRESSURE  EVOLUTION 
DURING  DECOMPOSITION 

In  Figure  1  we  show  the  time  evolution  of  the  potential  energy 
for  liquid  and  crystalline  nitromethane.  Figure  la,b,c  corresponds 
to  T  =  3000  K,  T  =  2500  K,  and  T  =  2000  K,  respectively.  After  a 
temperature  dependent  induction  time  during  which  the  energy 
remains  essentially  constant  and  the  first  chemical  reactions 
occur  (as  will  be  discussed  in  the  next  section)  the  potential 
energy  decreases  abruptly  as  exothermic  reactions  begin. 

In  all  cases  shown  in  Figure  1  the  liquid  samples  exhibit  faster 
chemistry  than  crystalline  ones;  however,  the  time  scales  of 
decomposition  in  our  simulations  are  not  well  converged  with 


respect  to  simulation  size.  Simulation  cells  large  enough  to  allow 
multiple  independent  initiation  events  are  necessary  to  obtain 
accurate  decomposition  time  scales.  We  are  currently  assessing 
the  cell  size  dependency  of  reaction  times. 

In  Figure  2  we  show  the  time  evolution  of  the  pressure  during 
the  reaction  of  the  liquid  1  sample  at  temperatures  T  =  2000, 
2500,  and  3000  K.  Before  the  chemical  reactions  occur  the 
pressures  range  from  ~18  GPa  at  T  =  2000  to  22  GPa  at  T  = 
3000  K.  As  the  reactions  proceed  and  the  number  of  molecules 
increases,  the  pressures  increase  by  ~  15%;  the  final  pressures 
range  from  ~23  to  ~29  GPa. 

In  order  to  validate  the  accuracy  of  the  ReaxFF  description 
of  mechanical  properties  under  such  large  compression  and 
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Figure  3.  Time  evolution  of  the  temperature  (a)  and  pressure  (b)  of 
crystalline  nitromethane  during  a  fast  annealing  from  T  =  1000  K  to  T  = 
4000  K.  QM-based  MD4  at  T  =  1000  K  lead  to  a  pressure  of  25  GPa  in 
good  agreement  with  our  results  (24.7  GPa).  When  the  temperature  is 
increased  to  4000  K  (in  a  ~0.6  ps  run)  the  pressure  rises  to  35  GPa,  also 
in  agreement  with  our  results.  Note  that  for  longer  times  our  simulations 
predict  a  further  increase  in  the  pressure. 


temperature  we  performed  a  simulation  at  a  higher  density 
(p  =  2.2  g/ cm3)  since  pressures  from  QM-based  MD  simulations 
are  available  for  this  density.4  As  in  ref  4  we  performed  a  short 
simulation  (0.2  ps)  at  T  =  1000  K  and  then  quickly  heated  the 
system  up  to  T  =  4000  K.  Using  ReaxFF  we  obtain  a  pressure 
P  =  24.70  GPa  for  the  T  =  1000  K  simulation  (averaged  over  the 
last  100  fs  of  the  simulation);  this  value  is  in  excellent  agreement 
with  the  reported  QM  result  (25  GPa).  In  Figure  3  we  show  the 
temperature  (a)  and  pressure  (b)  evolution  during  the  annealing 
run  that  takes  the  system  from  T  =  1000  K  to  T  =  4000  K.  The 
QM  results  indicate  that  the  pressure  rises  and  stabilizes  at  about 
35  GPa  (the  total  length  of  the  run  was  0.6  ps);  this  value  is 
indicated  as  a  horizontal  line  in  Figure  3b.  Again  the  ReaxFF 
results  are  in  very  good  agreement  with  the  QM  results.  Note  that 
for  times  longer  than  the  QM  run  the  pressure  increases  further 
to  about  40  GPa.  These  results  indicate  that  ReaxFF  provides  a 
very  accurate  description  of  the  equation  of  state  of  nitromethane 
at  high  pressures  and  temperatures. 

■  DECOMPOSITION  AND  REACTION  CHEMISTRY 

The  reactive  MD  simulations  in  this  paper  provide  a  very 
detailed,  molecular-level  description  of  the  decomposition  and 
reaction  of  nitromethane  under  various  conditions.  Such  infor¬ 
mation  should  allow  one  to  extract  valuable  information  about 
the  complex  chemistry  involved,  including  uni-  and  multimole- 
cular  reactions,  catalytic  steps,  etc. 

In  Table  1  we  show  the  times  corresponding  to  the  first  bond 
breaking  and  formation  for  a  variety  of  bonds  and  for  all  our 
simulations.  In  each  case  the  first  bond  breaking  event  is  marked 
with  an  asterisk.  Our  simulations  show  that  the  first  reaction 
during  decomposition  is  either  an  intermolecular  proton  transfer 
leading  to  CH2NO2  and  CH3NOOH,  or  an  isomerization 
reaction  leading  to  methyl  nitrite  (CH3ONO).  The  proton 
transfer  involves  breaking  a  C— H  bond  and  forming  a  new  bond 
between  the  proton  and  an  oxygen  in  a  nearby  molecule  (notice 
the  correlation  between  C— H  breaking  and  O— H  formation 
times  in  Table  l).  The  partial  charge  of  the  H  atom  transferred 
during  the  reaction  is  approximately  +0.2  e,  indicating  that  under 
such  condensed  phase  conditions  the  reaction  is  not  ideally  ionic. 


Table  1.  Times  (in  ps)  of  First  Bond  Breaking,  Bond  For¬ 
mation,  and  Intramolecular  Proton  Transfer  for  the  Various 
Simulations” 


bond 

bond  breaking  formation  intramolecular 


configuration 

C-H 

C-N 

N-O 

O— H  C— O  proton  transfer 

crystal  (T  = 

3000  K) 

0.1* 

2.8 

1.6 

0.1 

0.1 

1.38 

liquid  1  (T  = 

=  3000  K) 

0.05* 

0.5 

0.4 

0.05 

0.1 

1.13 

liquid  2  (T  = 

=  3000  K) 

0.3* 

0.7 

0.7 

0.3 

0.5 

1.17 

crystal  (T  = 

2500  K) 

2.4 

0.7* 

2.4 

2.4 

1.0 

liquid  1  (T  = 

=  2500  K) 

0.9 

0.6* 

1.5 

0.9 

0.6 

liquid  2  (T  = 

=  2500  K) 

0.5* 

2.2 

1.5 

0.4 

0.7 

0.45 

crystal  (T  = 

2000  K) 

28.0 

26.4* 

33.1 

28.0 

26.5 

51.95 

liquid  1  (T  = 

=  2000  K) 

2.0* 

7.4 

2.4 

2.0 

2.0 

37.42 

liquid  2  (T  = 

=  2000  K) 

4.6 

0.7* 

5.1 

4.6 

0.9 

a  Asterisks  mark  first  bond  breaking  events. 


M-  yk 

Figure  4.  Snapshots  of  isomerization  reaction  leading  to  methyl  nitrate 
obtained  from  the  MD  simulation  of  liquid  1  at  T  =  2000  K.  The  central 
snapshot  corresponds  roughly  to  the  transition  state  whose  structure  is 
in  good  agreement  with  recent  ab  initio  results.”8 


This  is  in  agreement  with  Mulliken  charges  in  QM-based 
simulations  of  Manaa  et  al.;4  they  find  a  charge  of  about  +0.3e 
for  the  same  atom.  In  Figure  4  we  show  five  snapshots  of 
the  isomerization  reaction  obtained  from  our  MD  simulation 
of  the  liquid  1  structure  at  T  =  2500  K;  the  time  elapsed  between 
the  second  and  the  fourth  snapshots  is  20  fs.  The  central  structure 
depicts  roughly  the  transition  state.  It  is  worth  emphasizing  that, 
while  a  single  molecule  is  shown  in  Figure  4,  our  simulation 
resolves  in  full  atomistic  detail  the  dynamics  of  32  molecules  in 
the  condensed  phase  at  high  density;  all  multimolecular  effects 
are  resolved  exactly  in  our  simulation.  The  C— N  and  C— O  bond 
distances  corresponding  to  the  condensed  phase  transition  state 
we  observe  with  ReaxFF  are  1.89  and  1.60  A,  respectively.  There 
has  been  some  controversy  in  the  literature  regarding  the 
structure  of  the  transition  state  between  nitromethane  and 
methyl  nitrate,  see  for  example,  ref  28.  Some  authors  have 
predicted  a  loose  transition  state  (with  C— N  and  C— O  bond 
distances  longer  than  3  A)29  while  others  predict  a  tight  transi¬ 
tion  state  with  bond  distances  around  2  A.  Our  MD  simulations 
predict  a  tight  transition  state  in  agreement  with  recent  ab  initio 
results;28  we  did  not  observe  an  initial  reaction  matching  the 
loose  transition  state.29 

Although  more  simulations  would  be  necessary  to  obtain 
better  statistics,  our  simulations  indicate  that  the  higher- 
frequency,  higher-activation  energy  C— H  breaking  process  is 
more  likely  to  be  the  initial  reaction  at  high  temperatures 
whereas  the  lower-frequency  lower-activation  energy  isomer¬ 
ization  reaction  leading  to  methyl  nitrate  tends  to  be  the  initial 
reaction  at  low  temperatures.  This  is  in  qualitative  agreement 
with  ref  29. 

All  the  first  chemical  events  shown  in  Table  1  happen  during 
the  induction  time  where  the  potential  energy  is  approximately 
constant;  triangles  in  Figure  1  mark  the  time  of  the  first  N— O 
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Figure  5.  Time  evolution  of  nitromethane  (green  lines),  several  important  products  [H20  (black),  C02  (red),  N2  (magenta),  and  NH3  (blue)],  and  the 
intermediate  H3CNOOH  (dashed  lines)  for  the  liquid  and  crystalline  structures  at  T  =  3000  K  and  T  =  2000  K.  Populations  are  given  per  nitromethane 
molecule. 


bond  breaking  process.  We  see  that  the  exothermic  chemistry 
starts  not  long  after  the  first  N— O  bond  breaking. 

Using  QM-based  MD,  Manaa  and  collaborators4  studied  the 
decomposition  of  solid  nitromethane  heated  to  3000  K  and 
found  that  the  first  reaction  was  an  intermolecular  proton 
transfer.  The  resulting  time  scales  for  this  first  reaction  were 
0.059,  0.780,  and  0.785  ps  in  three  independent  simulations. 
Our  results  (0.05,  0.1,  and  0.3  ps  also  in  three  independent 
simulations)  are  in  quantitative  agreement  with  the  QM  results; 
this  provides  an  important  direct  validation  of  ReaxFF  against  a 
first  principles,  QM-based  description  of  the  atomic  interactions. 

In  many  simulations  we  observe  intramolecular  proton  trans¬ 
fer  from  the  methyl  to  the  nitro  group  leading  to  the  aci  acid  form 
of  nitromethane  CH2NOOH;  Table  1  shows  the  correspond¬ 
ing  times  for  this  isomerization  reaction.  This  result  is  also  in 
agreement  with  QM  results.4 

For  the  crystal  at  T  =  3000  K,  we  find  that  the  very  first  water 
produced  contains  a  hydroxyl  group  from  the  nitrogen  end  of  a 
COH— NOOH  and  a  proton  from  a  formaldehyde.  The  hydroxyl 
group  is  formed  by  intramolecular  proton  transfer  from  the  CH3 
end  of  a  nitromethane  to  the  N02  end  of  another  one.  The 
carbon  end  of  that  molecule  subsequently  undergoes  several 
reactions,  resulting  in  COH— NOOH.  The  formaldehyde  is 
created  through  several  steps  beginning  with  the  dissociation 
of  a  CHj  group  from  a  nitromethane.  Subsequently,  this  group 
lost  a  hydrogen  and  then  gained  an  oxygen  by  inducing  a  N— O 
bond  to  break  in  a  nitromethane  remnant.  The  resulting  for¬ 
maldehyde  then  contributed  one  of  its  protons  to  form  water. 

In  Figure  5  we  show  the  time  evolution  of  the  population  of 
nitromethane  as  well  as  several  important  products  (H20,  C02, 
N2,  and  NH3)  and  the  intermediate  (H3CNOOH)  for  both 


Si-ping  Han  Fig  6 


Figure  6.  Bond  breaking  and  forming  rates  for  the  liquid  1  structure  at 
T  =  2000  and  T  =  3000  K. 


liquid  and  crystalline  structures  and  two  temperatures.  As  already 
mentioned  the  formation  of  CH3NOOH  via  a  proton  transfer  is  the 
first  reaction  that  appears  in  the  populations  plot.  The  formation  of 
H20  (the  main  product)  occurs  after  most  of  the  original  nitro¬ 
methane  molecules  have  undergone  the  initiation  reactions.  The 
final  population  of  H20  is  not  very  sensitive  to  temperature;  we  find 
about  0.8  H20  molecules  per  nitromethane.  Other  important 
products  are  N2,  C02,  and  NH3.  As  expected  the  population  of 
these  small  molecules  is  larger  at  T  =  3000  K  than  at  T  =  2000 IC;  in 
particular,  we  find  essentially  no  N2  molecules  at  T  =  2000  K. 

In  Figure  6  we  plot  the  rate  of  bond  breaking  and  formation 
(number  of  bonds  broken  or  formed  per  unit  time  per  nitro¬ 
methane  molecule)  for  various  key  pairs  of  atoms.  As  already 
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Figure  7.  Time  evolution  of  the  diffusion  constant  of  O  atoms  for  liquid 
and  crystalline  nitromethane  at  T  =  2000  K. 


mentioned  the  initial  reactions  involve  breaking  C— H  and  C— N 
bonds  and  forming  O  — H  and  C— O  ones.  We  can  see  from 
Figure  6  the  time  scales  of  the  various  reactions.  For  long  times 
we  find  the  main  events  are  O— H  bond  breaking  and  forming 
(with  equal  rates).  This  is  because  H20  is  the  main  product  and 
the  molecules  at  this  high  temperatures  undergo  frequent 
chemical  reactions  under  dynamical  equilibrium. 

Besides  the  small  molecules  reported  above,  decomposition 
leads  to  several  larger  molecular  fragments.  These  molecules 
contain  between  ~70%  and  ~40%  of  the  total  carbon  and 
nitrogen  of  the  system,  but  we  find  no  sign  of  initial  carbon 
clustering  within  the  relatively  short  time  of  our  simulations;  the 
largest  number  of  C  atoms  observed  in  a  molecule  is  three.  This  is 
consistent  with  experimental  observations. J° 

■  ATOMIC  DIFFUSION  DURING  REACTIONS 

Diffusion  of  atomic  species  during  the  decomposition  and 
reaction  of  energetic  materials  is  an  important  phenomenon  that 
governs  molecular  intermixing  and  plays  an  important  role  in  the 
degree  of  completion  of  the  chemical  reactions  during  detona¬ 
tion  or  combustion  as  well  as  in  the  process  of  carbon  clustering. 
The  formation  of  condensed-phase  aggregates  composed  mostly 
of  C  atoms  is  most  important  since  it  governs  the  balance 
between  CO  and  C02  (see  for  example  ref  12)  which,  in  turn, 
can  dramatically  affect  the  exothermicity  of  the  HE  material. 

Atomic  diffusivity  is  not  easily  accessible  experimentally  under 
the  extreme  conditions  characteristic  of  HE  materials,  making  the 
estimates  from  MD  very  valuable.  The  diffusion  constant  (D)  can 
be  defined  from  the  time  evolution  of  the  atomic  mean  square 
displacement  [Ar2(f)] 

Ar2(f)  =  Tlim  ijT  df'i£  [(r(t'  +  t)  -  r(t')}2  =  6Dt  (l) 

where  N  is  the  number  of  atoms  involved. 

In  describing  diffusion  it  is  important  to  wait  sufficiently  long 
times  such  that  Ar2  becomes  proportional  to  time  (Fickian 
behavior).  This  is  not  possible  in  our  case  since  the  decomposi¬ 
tion  occurs  in  relatively  short  times.  Consequently,  we  define 
an  effective,  local  in  time,  diffusion  constant.  We  obtain  the 
effective  diffusion  constants  from  our  MD  simulations  by  aver¬ 
aging  over  8  ps. 

In  Figure  7  we  show  the  time  evolution  of  the  diffusion 
constant  of  oxygen  atoms  at  T  =  2000  K  both  for  liquid  and 
crystalline  samples.  Despite  the  large  fluctuations  of  the  data  we 


find  the  initial  oxygen  diffusivity  in  the  liquid  samples  is  larger 
than  that  in  the  solid  one.  This  is  consistent  with  our  observation 
of  faster  initial  decomposition  of  the  liquid  samples,  see  Figures  1 
and  5.  The  decomposition  and  subsequent  chemical  reactions 
lead  to  a  decrease  in  the  molecular  size  of  the  system  and 
consequently  an  increase  in  diffusivity:  our  simulations  show 
an  increase  of  approximately  50%.  Finally,  as  the  products  are 
formed  information  about  the  initial  structure  of  the  material 
(crystalline  or  liquid)  is  lost  and  both  classes  of  samples  exhibit 
similar  diffusivities. 

■  DISCUSSION:  ReaxFF  VALIDATION 

Our  results  for  the  crystalline  sample  at  T  =  3000  K  are  in  good 
agreement  with  recent  QM-based  MD  simulations4  not  only 
regarding  the  initial  chemical  reactions  of  the  decomposition  and 
the  associated  time  scales  but  also  in  the  equation  of  state  of  NM 
up  to  very  high  pressures.  The  pressures  obtained  from  our 
simulations  at  T  =  1000  and  4000  K  and  p  =  2.2  g/ cm3  (24.7  and 
~35  GPa)  match  the  reported  QM  results  (25  and  35  GPa). 
Note  that  while  ab  initio  data  on  nitromethane  has  been  used  to 
develop  ReaxFF  (chemical  reactions  and  vibrational  properties), 
none  of  the  QM  MD  results  in  ref  4  were  included  in  anyway  into 
the  ReaxFF  parametrization. 

Reference  4  used  density  functional  theory  with  ultrasoft 
pseudopotentials,  a  plane  waves  basis  set,  and  the  spin-polarized 
generalized  gradient  approximation  of  Perdew  and  Wang  to 
describe  exchange  and  correlation.  This  level  of  DFT  is  generally 
very  accurate  providing  an  excellent  validation  of  ReaxFF 
which  has  been  parametrized  against  similar  quality  DFT  data 
(no  pseudopotential  with  the  B3LYP  flavor  of  DFT  including  the 
Becke  GGA  exchange  functional  with  the  LYP  correlational 
functional  and  some  exact  exchange). 

These  results  indicate  that  ReaxFF  accurately  describes  com¬ 
plex  chemical  processes  in  CHNO  systems  including  multi- 
molecular  reactions,  ionic  reactions,  and  mechanical  properties 
at  extreme  conditions.  They  are  similar  to  results  from  another 
set  of  recently  reported  ReaxFF  simluations  on  thermal  decom¬ 
position  of  liquid  NM.31  However,  our  results  include,  in 
addition,  simulations  on  crystalline  systems  and  a  much  more 
detailed  analysis  of  the  chemistry.  This  detailed  and  direct 
validation  of  the  accuracy  of  the  ReaxFF  description  of  the 
interatomic  forces  is  important  since  it  indicates  that  we  can 
use  ReaxFF  for  much  larger  systems  and  longer  times  for  which 
the  QM  MD  is  currently  impractical.  This  also  suggests  that  the 
functional  forms  in  ReaxFF,  including  the  description  of  charge 
transfer  and  of  bond  orders,  capture  the  physics  responsible  for 
the  chemistry.  In  addition  it  supports  the  value  of  using  an 
extensive  training  set  containing  over  1515  separate  reactions 
and  equilibrium  molecules  designed  to  characterize  atomic 
interactions  under  various  environments  [likely  and  unlikely 
(high  energy)],  including  over  50  molecules  describing  nitro¬ 
methane  and  its  reactions  with  water  and  NH3.  The  training  set 
contains  bond  dissociation  and  compression  curves  of  all  possi¬ 
ble  bonds,  angles,  and  torsion  bending  for  all  cases,  vibrational 
properties,  and  equation  of  state  data  of  crystals.j2 

■  CONCLUSIONS 

We  used  ReaxFF  with  MD  to  study  the  thermal  decomposi¬ 
tion  of  crystalline  and  liquid  nitromethane  at  high  temperatures 
(2000—3000  K)  and  high  density  (p  =  1.97  g/cm3).  At  the 
highest  temperature  (T  =  3000  K)  the  first  chemical  reaction  is 
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an  intermolecular  proton  transfer  process  leading  to  CH2N02 
and  CH3NOOH.  At  lower  temperatures  (T  =  2000  K  and  T  = 
2500  K)  the  first  reaction  is  often  an  isomerization  leading  to 
CH3ONO.  We  also  find  at  very  short  times  an  internal  proton 
transfer  isomerization  reaction  leading  to  aci-nitromethane 
(CH2NOOH).  Crystalline  and  liquid  samples  exhibit  the  same 
overall  chemistry.  The  main  decomposition  product  is  H2Oj  we 
find  ~0.8  H20  molecules  per  original  nitromethane  in  all  cases 
studied.  Other  products  we  observe  are  C02,  NH3)  and  N2. 
Finally,  while  the  time  scales  we  obtain  from  the  ReaxFF 
simulations  are  similar  to  previously  published  rates  in  simula¬ 
tions  using  quantum-mechanics-based  simulations,  the  reaction 
zone  in  NM  involves  longer  time  scales  of  a  few  nanoseconds,  see 
ref  30  for  a  compilation  of  data.  We  believe  an  important 
contributing  factor  for  this  discrepancy  is  that  during  detonation 
the  reactive  material  expands  as  the  reaction  proceeds  in  the 
reaction  zone  leading  to  slower  chemistry. 
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